Monte Carlo study of the elastic interaction in heteropitaxial growth 
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We have studied the island size distribution and spatial correlation function of an island growth 
model under the effect of an elastic interaction of the form 1/r 3 . The mass distribution P„(t) that 
was obtained presents a pronounced peak that widens with the increase of the total coverage of the 
system, 9. The presence of this peak is an indication of the self-organization of the system, since it 
demonstrates that some sizes are more frequent than others. We have treated exactly the energy of 
the system using periodic boundary conditions which were used in the Monte-Carlo simulations. A 
discussion about the effect of different factors is presented. 
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I. INTRODUCTION 

Epitaxial growth has been the focus of much interest 
in the past years. This interest is derived mainly from 
the fact that this kind of processes have numerous appli- 
cations in developing new types of devices and ma- 
terials with some special characteristics. Among these 
processes, there is one type of growth characterized by 
the presence of long-range elastic interactions that play 
a special role 0. Molecular Beam Epitaxy (MBE) un- 
der the effect of elastic interactions has been a topic of 
recent research, with many applications, specially in the 
fabrication of low dimensional nanostructures, like quan- 
tum dots (QD) and quantum wires 0^,0] ■ One chal- 
lenging question concerning this topic is how, and by 
which mechanism, island organization occurs. A 
wide range of material/substrate combinations have been 
observed (e.g., InAS on GaAs). When one kind of ma- 
terial is deposited over a different one -, usually, with a 
different lattice parameter - it will induce, through this 
structural difference, a long-range elastic interaction be- 
tween the deposited atoms during surface growth. The 
deformation thus obtained originates a strain in the sub- 
strate that causes different particles to repel each other. 
This kind of interaction is supposed to be the mechanism 
responsible of the self-organization observed experimen- 
tally. Many recent studies were developed in order to 
identify the influence of strain on epitaxial and surface 
morphology during growth j|,[ll]-[l6| . 

Several authors have studied the effect that strain in- 
duces on epitaxial growth on different types of systems. 
For instance, the effect of elastic strain on the proper- 
ties of the well know Eden model [|l6| , and for other ver- 
sions of a harmonic interaction between the lattice atoms 
jl7],[l8]]. Also, the Lennard- Jones potential was used to 
study a similar phenomena Jl9[ |. 

In this paper we shall consider that the strain induced 
on the system is due to an repulsive elastic interaction be- 
tween the deposited particles proportional to 1/r 3 . This 



type of potential can be derived from elasticity theory 
considerations ppc|-|2^|, when a lattice distortion is cre- 
ated (e.g. by cutting out a sphere of the bulk and sub- 
stituting it by a different radius sphere) a field of lattice 
strains is created. It is already known that this type of 
long-range interaction can be applied to the absorption of 
atoms onto a surface but is only valid on the case of very 
thin absorbed clusters (submonolayer regime); in more 
general cases it can be obtained by scaling laws (lq|. 



II. ELASTIC INTERACTION POTENTIAL 

We define the elastic potential to be of the form, 
E = GmiiTij /r 3 , where r represents the distance between 
the particles, mi is the "mass" of the particle i, and G is 
the coupling constant. G usually depends on the elastic 
properties of the substrate, such as the Young modulus, 
the Poisson ratio and the lattice mismatch. The coupling 
constant is given by Q as being, G = 7r(l — a 2 )a 2 f 2 / E, 
but other authors define it in different ways ^|{n|. The 
numerical values of the various forms differ on several or- 
ders of magnitude. We overcome this difficulty by leaving 
the analythical form of G unspecified and determining a 
reasonable numerical value for it by physical considera- 
tions. 

Since we are trying to study the behavior of a macro- 
scopic material, we employed periodic boundary condi- 
tions during the simulations. The presence of this type of 
boundary conditions implies that some treatment must 
be given to the energy defined above in order to avoid 
some undesirable "finite-size effects" that could originate 
"unphysical" results. To accomplish this, we consider an 
infinite succession of replicas of the system and calculate 
the total energy of the infinite system thus obtained. The 
total energy has the contribution of two components, the 
first being the interaction energy between the particle 
that is currently suffering the absorption on the origi- 
nal system and all the other copies of this particle that 
belong to the other systems. This contribution is given 
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where L is the linear dimension of the system and is 
the Riemman Zeta function. The first sum is performed 
over all particles present in the system at this time. 

The second contribution to the total energy is given 
by the interaction between the deposited particle and all 
particles deposited previously in the system. This con- 
tribution is expressed as, 
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where the first sum is performed over all pairs of particles, 
dij is the distance that separates the two particles in the 
original system divided by L. Considering that 
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where t])W(ol) is the second order Polygamma function, 
and after performing some straightforward algebraic ma- 
nipulations, the energy Ei takes the following form, 
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The total energy of the system is then given by 
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During the simulations, and since we are not interested 
in the absolute value of the energy but only in energy 
differences, we shall only consider the "effective" value 
of the energy, that is, the part of the energy that is not 
constant, in order to shorten the CPU time without loss 
of precision in the results. 



to the number of nearest neighbor (NN) sites that are 
occupied. When only one NN is occupied, the particle 
adheres irreversibly to the preexisting cluster and an- 
other deposition attempt is performed. If the two NN 
are occupied, the particle adheres to both clusters, co- 
alescing them to become one single cluster with mass 
conservation. Finally, if none of the NN sites is occu- 
pied, the particle diffuses, due to the repulsive effect of 
the potential generated by the mass distribution present 
in the system, moving away from the larger cluster and 
becoming closer to the smaller one, until it reaches the 
local minimum of the energy. At each diffusion step, the 
energy resulting from the interaction of the adatom with 
every other particle present in the system is calculated. 
At this point, the particle begins to diffuse due to the ef- 
fect of the temperature with probability proportional to 
e -AE/k B T ^ w h ere jg the total energy of the system 
During this process, a number D of random steps is per- 
formed. If during this random walk motion, the particle 
collides with another particle, it aggregates irreversibly 
and another particle is deposited. This model has sev- 
eral adjustable parameters such as the temperature T, 
the number of diffusion steps D and the value of the in- 
teraction constant G. In order to make the simulations 
behave realistically, these parameters must be adjusted 
and their effect on the final result must be studied and 
well understood. This has been done by varying the pa- 
rameters in order to see the effect that each parameter 
individually has on the final result. During the experi- 
mental study of this type of processes, one usually uses 
temperatures in the interval 300A < T < I000A'. To 
use this values during the simulations, one must adjust 
G in a way as to make the factor AE/ksT ~ 1. The 
typical value of the energy differences in this model is of 
the order 10 9 , and considering that T ~ 10 2 we find that 
Gm 2 AE ef /(L 3 k B T) ~ 1, and so, Gm 2 /L 3 k B ~ 1(T 7 . 



IV. RESULTS 

In this section we will present the results obtained us- 
ing Monte Carlo simulations. To characterize the coars- 
ening dynamics, two quantities were sampled and aver- 
aged over the initial conditions: cluster "mass" distribu- 
tion, P n (t) and the correlation function at equal times 
C(r, t), defined by 



III. SIMULATIONS 



C(r,t) = {p(r' + r,t)p(r',t))~{p(r',t)) 2 



(0) 



The model described earlier was implemented in a rel- 
atively simple way. We consider the substrate to be one 
dimensional and we shall only consider the regime of sub- 
monolayer growth. One site of the system is selected ran- 
domly. If that site is occupied, the deposition attempt 
fails and another site is selected. If the selected position 
is empty, three possible situations can occur according 



where p(r, t) is the site density. The most convincing re- 
sult yielding the self-organization process is the fact that 
the "mass" distribution function presents a well defined 
peak. 
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FIG. 1. Mass distribution variation with the increase of the 
system coverage. Notice that the distribution function widens 
as the coverage increases. In the inset we represent the num- 
ber of clusters present in the system as a function of coverage. 
It is visible that after a certain point, the number of clusters 
decreases due to the occurrence of coalescence phenomena. In 
all the simulations, T = 500K and D = 200 steps. 



The shape of the distribution is maintained as one in- 
creases the coverage, but the height of the function tends 
to decrease as the width increases. This fact was ex- 
pected to happen, larger values of coverage imply that 
fewer individual clusters are present in the system, but 
with large sizes. 

There exist two basic process that the system has avail- 
able in order to organize itself. The first, the nucleation 
of new clusters, is dominant in the early stages of the 
system evolution, when the coverage is small and the 
adatoms never collide. The second, is the coalescence 
of existing clusters. This process becomes dominant as 
the coverage increases, originating larger clusters but in a 
smaller number. This two regimes are clearly seen in the 
inset of Fig.[j]. In the beginning, the number of clusters in 
the system seems to grow almost linearly with the cover- 
age. Afterwards, there exists a crossover period when the 
number of clusters is approximately constant, that hap- 
pens when the growth of existing clusters becomes more 
frequent than the nucleation of new ones. Finally, coales- 
cence begins to dominate the dynamics and the number 
of clusters in the system diminishes until it becomes one 
when the coverage gets very large. 

The behaviour of the system is described by the island- 
size distribution function, P n (9). Assuming there exists 
a scaling for P n {9) one may write, P n (6) = A/(A Q n, \°9). 
The coverage 9 grows with time, but satisfies, 9 = 
Sn>i n Pn{&)- This sum can be approximated by an in- 
tegral, 9 w J °° nP n (9)dn resulting a relation between the 
exponents, = 2a — 1. Therefore, one can write, 
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All this results follows from the assumption that there ex- 
ist only one characteristic size in the system, the average 



island size, S = £ n >! nP n (0)/E„>i Pn(0) ~ a/P - The 
data colapse of P n (t) is shown in Fig.||. It was obtained 
for a = 1.47 ± 0.05 and (3 = 1.23 ± 0.05. 
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FIG. 2. Scaling of the mass distribution function Pn(t) 
with the system coverage. The values of a and j3 are re- 
spectively, 1.23 and 1.47 with an error of 0.05 in both cases. 
(T = 500A' and D = 200) 
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FIG. 3. Variations in the mass distribution when the num- 
ber of diffusion steps D and the temperature T changes, a) 
With the increase in D, the mass distribution keeps the same 
basic shape, but it's width increases while it's height decreases 
in a form similar to the one observed in Fig^] when we in- 
creased the coverage of the system (T = 500 and 9 = 0.3). 
The temperature doesn't cause any important changes in the 
form that the system organizes itself.(T = 500K, 9 = 0.3) 

When we keep the value of the coverage fixed and vary 
D, the mass distribution function behaves in a manner 
similar to the one described above. As D increases, the 
distribution function widens and flattens. This is due to 
the fact that, with a larger number of diffusion steps, the 
adatom has a larger probability of diffusing away from 
the local minimum of energy and coalescing with other 
particles present in the system, originating larger clus- 
ters. On the other hand, when we change the tempera- 
ture, nothing seems to happen, the distribution function 
maintains its basic properties. 
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FIG. 4. Correlation function at various coverages. In a) 
we observe the evolution of the correlation function as the 
coverage increases. It is clearly seen that the position of the 
minimum tends to move to larger values of r (T = 500K 
and D = 200). b) A similar precession of the minimum is 
observed when D is increased (9 = 0.3 and T = 500). c) 
The correlation function doesn't seem to be affected by the 
changes in T. This is probably due to the fact that only the 
diffusing adatom feels the effect of the temperature (D = 200 
and 9 = 0.3). 

As it can be seen in Figfil a, the correlation function 



C (r, t) , defined in Eq|| displays a characteristic behav- 
ior, starting at values of the order of 0.2 and decreasing 
until a minimum, at distances of the order of ten lattice 
units, and finally oscillating with decreasing amplitude 
around zero until it becomes effectively zero at distances 
of the order of 40. As it is easily seen from FigJ^, the cor- 
relation function always maintains the same shape, even 
when 9 or D are increased. In the latter case, the posi- 
tion of the minimum seems to move to larger distances 
as D is increased, which means that the system becomes 
more correlated. Once again, the temperature doesn't 
have any real effect on the results. 

A parameter that usually has a great importance in 
experimental study of this type of systems is the tem- 
perature. As can be seen in the previous figures, the 
temperature doesn't seem to have a great influence in 
the final result, contrary to what was expected. This pe- 
culiar behavior of the system can probably be explained 
by the fact that only the adatom that is being currently 
absorbed feels the temperature, the rest of the system 
being actually frozen. If we allowed the system to re- 
arrange itself after the deposition of each particle, the 
temperature dependence would probably be more real- 
istic, but the computation time required would also be 
much higher. 



In conclusion, the simulations carried in a one dimen- 
sional system in the submonolayer regime with long range 
interactions allowed us to observe the mechanism of self- 
organization through the formation of islands of similar 
size over all system. The influence of different factors 
on this behavior was tested. Surprisingly, it does not 
show any dependence on the temperature, at least for 
the tested range (300 A" < T < 1000 A). 

Ordering occurs to minimize the repulsive elastic inter- 
actions between absorbed atoms. This self-organization 
breaks down when the coverage gets large which makes 
the adatom have less space to find an equilibrium posi- 
tion and makes the coalescence events become more and 
more frequent and finally dominate the dynamics of the 
system. We are now extending this results to the 2D 
case. 
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